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ABSTRACT 

An  implicit  pressure  correction  method  on  unstructured  Cartesian  grid  is  developed  for  the 
incompressible  Navier-Stokes  equations.  An  immersed  boundary  method  is  also  incorporated  to  treat 
the  body  geometry.  The  pressure  Poisson  equation  is  solved  by  a  multi-grid  method.  A  fourth  order 
artificial  dissipation  is  added  to  the  pressure  field  to  suppress  the  even-odd  decoupling.  Tests  show 
that  with  an  appropriate  amount  of  dissipation,  the  method  is  second  order  accurate  both  in  time  and 
space.  The  driven  cavity  flows  with  and  without  immersed  bodies  are  computed  to  demonstrate  the 
capability  of  the  present  scheme. 

Keyword:  Implicit  Pressure  Correction  Method,  Unstructured  Cartesian  Grids,  Immersed  Boundary 
Method 

1.  Introductions 

The  numerical  method  for  the  incompressible  Navier-Stokes  equations  has  been  highly  developed 
in  recent  years.  One  popular  method  is  the  pressure  correction  method  (or  fractional  step  method)  [1-5] 
in  which  the  divergence-free  condition  is  obtained  via  an  appropriate  pressure  field.  In  these  methods 
the  velocity  is  updated  by  a  2-step  predictor-corrector  time  integration.  In  the  predictor  step  the 
intermediate  velocity  field  computed  from  the  momentum  equations  is  usually  not  divergence  free. 
Then  in  the  corrector  step  a  pressure  field  is  applied  to  correct  the  non-divergence-free  part  of  the 
intermediate  velocity.  This  pressure  field  obeys  a  Poisson  equation.  Normally  in  these  methods 
staggered  grids  are  used  to  achieve  pressure-velocity  coupling  [1-3].  If  a  non-staggered  grid  is  used, 
then  a  cell-face  velocity  is  defined  in  addition  to  the  cell-center  variables  [4,5]  to  compute  the  volume 
flux.  The  time  integration  method  used  for  velocity  is  usually  semi- implicit.  That  is,  the  predictor 
step  of  the  time  integration  usually  employs  an  explicit  scheme  (Adams-Bashforth  or  RK3)  for  the 
inviscid  terms  and  an  implicit  scheme  (Trapezoidal  rule)  for  the  viscous  terms.  It  is  well  known  that 
the  RK3  scheme  has  a  CFL  limit  of  -J~3  and  the  Adams-Bashforth  scheme  is  weakly  unstable  for  the 
model  convection  equation  [1].  These  stability  limitations  may  become  sever  when  the  grid  is  locally 
refined  or  stretched. 

Traditionally  a  body-fitted  [1,5]  structured  or  unstructured  grid  system  is  employed  to  handle  the 
immersed  bodies  in  the  flow  field.  More  recently  the  immersed  boundary  method  [2,3,4]  on  fixed 
Eulerian  grids  has  gained  popularity  for  their  ease  in  treating  the  complex  geometries.  In  Fadlun  et  al. 
[2]  and  Kim  et  al.  [3],  a  body  in  the  flow  field  is  modeled  by  adding  momentum  forcing  to  the 
appropriate  grid  cells  to  simulate  the  appropriate  boundary  condition.  In  Mittal  et  al.  [4],  the  grid  cells 
cut  by  the  body  are  treated  with  special  differencing  schemes.  The  usefulness  of  these  methods  lies  on 
the  appropriateness  and  simplicity  of  the  procedure  to  model  the  true  body  geometry  and  its  effects  on 
the  flow  field.  In  this  aspect,  the  capability  of  local  refinement  around  body  boundaries  is  desirable. 

In  this  work  an  implicit  pressure  correction  method  is  developed  to  solve  the  incompressible 
Navier-Stokes  equations.  The  fully  implicit  time  integration  is  achieved  by  the  technique  of  sub¬ 
iteration  [6].  The  unstructured  Cartesian  grid  system  is  employed  to  allow  local  refinement  when 
necessary.  The  pressure  Poisson  equation  is  solved  by  an  implicit  multi-grid  method  [7].  A  fourth 
order  artificial  dissipation  is  added  to  the  pressure  field  to  suppress  the  even-odd  decoupling  of 
pressure.  An  immersed  boundary  method  is  also  developed  to  handle  the  immersed  bodies.  The  no- 
slip  boundary  condition  is  enforced  by  the  method  of  direct  forcing  [2].  The  driven  cavity  flows  with 
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immersed  bodies  are  computed  to  demonstrate  the  capability  of  the  present  scheme. 


2.  Implicit  Pressure  Correction  Method 

The  integral  form  of  incompressible  Navier-Stoked  equations  can  be  written  as 
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where  v  and  P  are  Cartesian  velocity  and  pressure,  Re  is  Reynolds  number,  fB  is  the  momentum 
forcing  simulating  the  solid  body,  V  is  the  volume  of  the  control  volume  CV,  S  is  the  surface  CS  of 
the  cell,  R  represents  the  surface  integral  of  inviscid  and  viscous  fluxes  except  the  pressure  term. 
Assume  that  the  flow  state  at  time  level  n  is  known  on  a  cell-centered  non-stagger  grid.  An  implicit 
time  integration  with  sub-iteration  [6]  applied  to  the  momentum  equation  can  be  written  as: 
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where  the  superscript  n  is  the  time  level  index,  s  is  the  sub-iteration  index  for  velocity,  m  is  the  sub¬ 
iteration  index  for  pressure,  At  is  the  time  increment,  AV  is  the  cell  volume,  AS  is  the  cell  face,  the 
summation  operator  represents  a  closed  surface  integral  over  the  cell  faces,  and  0  <0  <1.  Equation  (2) 
is  iterated  in  s  with  fixed  m  and  n.  The  initial  conditions  are  Pm  =  P"  at  m=l  and  v  '  —vn  at  s=L 
When  the  sub-iteration  in  s  converges,  we  obtain  an  intermediate  state  v*  satisfying 
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Generally  this  v‘"  does  not  satisfy  the  divergence -free  condition.  A  pressure  correction  is  sought  to 
obtain  a  divergence-free  velocity  vm+I  and  pressure  Pm+I  that  satisfy  the  following  equations: 
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Div(ym+I)  =  0. 

where  Div(  •  )  represents  an  appropriate  divergence  operator.  Set  the  pressure  correction  as 
(J)  =  p"’+I  -  p"'  and  subtract  Eq.  (3)  from  Eq.  (4),  we  obtain  the  equation  for  velocity  correction: 
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Take  the  divergence  of  Eq.  (5),  we  obtain  the  Poisson  equation  for  the  pressure  correction: 
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Equation  (6)  is  to  be  solved  for  (j>  and  the  pressure  P"1+I ,  which  in  turn  is  to  be  substituted  into  Eq.  (5) 


to  obtain  the  divergence-free  velocity  vm+1 .  Equations  (2),  (5)  and  (6)  constitute  one  pressure 
iteration  in  m.  This  iteration  in  m  is  considered  converged  when  the  pressure  and  velocity  corrections 
in  Eq.  (5)  are  small  enough.  When  this  is  achieved,  the  latest  computed  flow  field  is  the  flow  state  at 
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the  new  time  level  n+1,  that  is,  v"+I  =  v"'+/  =  v  and  Pll+I  =  P",+I  =  Pm .  In  this  work,  Eq.  (2)  is 
inverted  by  an  approximate  LU  factorization  scheme  [6]. 


3.  Poisson  Equation  and  Cell-Face  Velocity 

To  evaluate  the  fluxes,  the  variable  values  at  the  cell  faces  need  to  be  reconstructed.  First,  the 
cell-vertex  values  are  obtained  by  a  distance-weighted  averaging  of  the  surrounding  cell-center  values. 
The  averaging  constant  is  proportional  to  the  inverse  of  the  distance  from  the  cell  center  to  the 
particular  vertex.  The  variable  gradients  can  be  estimated  by  differencing  the  cell-vertex  values 
enclosing  the  particular  cell.  Then  the  left  and  right  states  at  a  cell  face  can  be  linearly  interpolated 
using  the  gradient  just  obtained.  The  mean  cell-face  velocity  and  pressure  are  obtained  by  a  simple 
averaging  of  the  left  and  right  states: 

Vf,mean  =  0-5(v!f  +  vRf ),  PLmean  =  0.5{P^  +  Pf  )  (7) 

where  the  superscript  L  and  R  respectively  indicate  the  left  and  right  states  of  a  face,  the  subscript  / 
indicates  a  cell  face,  the  subscript  mean  indicates  the  mean  values.  The  mean  cell-face  pressure  is 
used  to  compute  the  surface  integral  of  pressure  in  Eq.  (2).  This  is  equivalent  to  a  central  differencing 
for  the  pressure  gradient  on  a  Cartesian  mesh.  The  mean  cell-face  velocity  v*fmean  is  used  to  compute 


the  divergence  of  v*in  Eq.  (6).  This  is  equivalent  to  a  central  differencing  for  the  divergence  operator. 

The  left  hand  side  operator  of  Eq.  (6)  is  equivalent  to  a  discretized  Laplacian.  However,  it 
becomes  non-compact  when  the  divergence  operator  is  applied,  making  the  inversion  of  Eq.  (6) 
expensive.  Thus,  for  computational  efficiency  the  left  hand  side  of  Eq.  (6)  is  simplified  as 
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The  left  hand  side  operator  is  a  compact  Laplacian  discretization  when  the  cell-face  V(J>  is  obtained  by 


differencing  the  neighboring  cell-center  values.  In  this  work,  Eq.  (8)  is  solved  by  a  multi- grid  method 
on  unstructured  Cartesian  grid  [7]. 

Note  that  Eq.  (5)  is  an  equation  for  cell-center  velocities  while  the  pressure  at  cell  faces  are 
interpolated.  But  Eq.  (8)  is  an  equation  for  cell-centered  (f)  with  v*f  mean  interpolated.  Because  of  this 


inconsistency  in  variable  locations,  the  cell-center  v"'+/ computed  by  Eq.  (5)  may  not  satisfy  the 
divergence- free  condition  exactly.  We  resolve  this  inconsistency  in  a  way  similar  to  the  approach  in 
Ref.  [4,5],  by  defining  a  new  cell-face  velocity.  That  is,  while  Eq.  (5)  is  still  used  to  update  the  cell- 
center  velocity  v",+1 ,  a  cell-face  velocity  v™+/  is  updated  from  mean  by 


-L - A^  +  0V4,  =  O  (9) 

At 

where  V(j>  at  the  cell  face  /is  obtained  by  differencing  the  neighbouring  cell-center  values.  It  is  easy 
to  verify  that  Eq.  (9)  is  consistent  with  Eq.  (8)  and  consequently  v'"+I  satisfies  the  divergence-free 


condition.  This  v"l+I  is  used  to  compute  the  volume  flux  at  the  next  iteration  in  m.  The  face-center 
v'”+1  is  updated  in  parallel  to  the  cell-centered  v"'+1 ,  while  the  flow  states  are  still  the  cell-center  v"'+1 
and  Pm+I . 


3.  Convective  and  Viscous  Fluxes 

When  iterating  Eq.  (2)  in  s,  the  last  past  divergence-free  velocity  V"  is  used  to  compute  the 
volume  flux.  Thus,  the  convective  momentum  fluxes  in  Rs  of  Eq.  (2)  are  computed  as: 

Vf-A S  (10) 

where  the  subscripts  L/R  represents  the  velocity  upwinding: 
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Equation  (10)  is  equivalent  to  a  second-order  accurate  upwind  differencing  for  the  convective 
momentum  fluxes.  The  viscous  fluxes  in  Rs  of  Eq.  (2)  are  computed  as 
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The  gradient  of  velocity  components  at  the  cell  faces  can  be  obtained  by  differencing  the  neighbouring 
cell-center  values.  On  a  regular  Cartesian  grid  this  is  equivalent  to  a  compact  second-order  accurate 
central  differencing  for  the  viscous  terms. 

5.  Artificial  Dissipation 

Our  experience  indicates  that  the  velocity  upwinding  and  the  use  of  v"‘  to  compute  the 

convective  fluxes  are  not  enough  to  prevent  the  even-odd  decoupling  of  the  pressure  field.  Thus,  an 
artificial  dissipation  is  added  explicitly  for  stabilization.  The  dissipation  for  a  cell  is  defined  as 

7disp  _  cdisp 


paisp  _ 


1  +  A x! Re 


Y,(PfR~Pf)n-e 


(13) 


cs 


where  n  is  the  unit  normal  of  a  cell  face  pointing  outward  of  the  cell,  e  is  the  Cartesian  unit  vector 
parallel  to  n  ,  and  we  take  0.01  <  cdisp  <0.1  in  this  paper.  This  dissipation  is  added  explicitly  to  the 
pressure  field  after  solving  Eq.  (8): 

(14) 
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It  is  easy  to  show  that  the  dissipation  is  proportional  to  Ax 


dx' 


for  a  Cartesian  mesh.  Our 


experience  indicates  that  with  an  appropriate  amount  of  cdisp,  the  added  dissipation  is  beneficial  and 
sometimes  crucial  to  obtaining  a  smooth  pressure  distribution. 


6.  Order  Analyses 

All  examples  shown  in  this  paper  are  computed  on  a  personal  computer  using  32-bit  single 
precision.  The  analytical  solutions  of  decaying  vortices  [8]  are  used  for  order  analyses: 

u(x,y,t )  =  -cos(a7tx)sin(a7tv)e  ~a  n  "Re ,  v(x, y,  t)  =  sin(a7tx)cos(a7tv)e  71  t,Re 


p  ^  _  cos (2am)  +  cos(2<my)  ^4a^t/Re 

tx,y,  J  4  e 

Here  we  take  -  8  <  x,y  <8  ,  a=l/8,  and  Re=l.  Exact  solutions  are  used  as  initial  and  boundary 
conditions.  The  Trapezoidal  rule  with  0  =0.5  is  used  for  time  integration.  The  iteration  in  m  is 
considered  adequate  when  both  the  L2  norms  of  divergence  field  and  the  pressure  correction  (f)  are 


less  than  7.x  10~5  and  at  least  one  of  the  two  is  less  than  0.5 x  / (T1  .  The  allowable  number  of 
pressure  iteration  is  set  to  5  <m<  30 . 

To  test  the  spatial  accuracy  of  the  scheme,  the  system  is  integrated  48  time  steps  to  t=0.3  on 
Cartesian  meshes  of  cell  length  Ax  =0.125,  0.25  and  0.5.  The  time  step  is  A t  =  0.00625  to  reduce  the 
time  discretization  error.  The  constant  cdisp=0.03  is  used.  The  L2  and  maximum  norms  of  the  error 
in  pressure  and  velocity  are  plotted  in  Fig.  1 .  The  average  order  of  accuracy  estimated  for  pressure  is 
2.32  in  L2  norms  and  1.86  in  maximum  norms.  For  velocity,  the  accuracy  is  2.03  in  L2  norms  and 
2.27  in  maximum  norms.  For  comparison,  similar  curves  are  plotted  for  cdisp=0  in  Fig.  2.  Without 
artificial  dissipation  in  pressure,  only  the  L2  norm  of  velocity  has  an  average  order  of  accuracy  of  1.88, 
all  other  norms  deteriorate  to  first  order  or  even  worse.  This  has  demonstrated  the  need  for  artificial 
dissipation  in  pressure. 

To  test  the  time  accuracy,  the  system  is  integrated  to  t=0.8  on  a  Cartesian  grid  with  Ax  =0.125. 
The  time  steps  are  At  =0.2,  0.4  and  0.8.  The  maximum  allowable  number  of  pressure  correction  h  m 
for  each  time  step  is  increased  to  1000.  The  velocity  and  pressure  norms  verse  time  steps  are  plotted 
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in  Fig.  3.  The  average  order  of  accuracy  for  all  norms  in  Fig.  3  is  second  order  or  higher.  The 
maximum  CFL  number  observed  in  the  calculations  for  Fig.  3  is  6.4.  This  indicates  that  the  sub¬ 
iteration  technique  has  eased  the  limitation  on  the  time  step. 

7.  Driven  Cavity  Flows 

The  numerical  solutions  of  a  driven  cavity  obtained  by  Ghia  et  al.  [9]  are  used  to  test  the  steady-state 
computation  of  the  present  method.  A  Cartesian  grid  of  size  128x128  is  used  to  discretized  a  square 
cavity  of  unit  length  in  size.  The  upper  wall  is  moving  to  the  right  at  unit  speed.  The  Euler  implicit 
method  with  0=7  is  used  with  m  <1  at  each  time  step.  We  take  cdisp=0.01 .  The  Reynolds  number 
Re  is  based  on  wall  speed  and  the  cavity  width.  For  the  case  of  Re=3200,  the  maximum  CFL  number 
is  set  to  2.  The  velocity  component  u(y)  along  the  vertical  centerline  and  the  velocity  component  v(x) 
along  the  horizontal  centerline  are  plotted  in  Fig.  4.  The  comparison  with  the  result  of  Ghia  et  al.  [9] 
who  also  used  a  similar  Cartesian  grid  is  very  good.  The  computed  path  lines  are  plotted  in  Fig.  5, 
showing  one  primary  vortex  at  the  center  and  three  secondary  vortices  on  the  corners  and  the  left  wall. 


8.  Immersed  Boundary  Method 

In  this  work,  we  use  a  level  set  function  ())/  defined  on  the  cell  vertices  to  track  the  presence  of 
the  immersed  bodies.  W  chose  to  be  a  signed  distance  function  whose  absolute  value  equals  the 


shortest  distance  to  the  body  surface.  It  is  defined  (| ),  >  0  outside  the  body,  (] );  <0  inside,  and  (j)7  =  0 
on  the  body  surface.  With  known,  the  volume  fraction  occupied  by  the  immersed  body  in  a  cell,  or 
the  Volume  of  Body  (VOB)  function  <f)6 ,  can  be  computed  easily.  We  have  <f)z,  =7  for  body  cells, 


§b  =0  for  fluid  cells  and  0  <§h  <  1  for  interface  cells  containing  the  body  contour  <\>,  =  0 .  We  take 
<f)ft  =  0.5  as  the  representative  body  surface.  The  accuracy  of  <\)h  depends  on  the  assumed  shape  of 
the  body  contour  of  ())/  =0  inside  the  interface  cells.  Here  we  assume  (| ),  =  0  contour  is  linear  inside  an 
interface  cell.  This  simplifies  the  calculation  of  ()/; ,  but  the  resulting  (j)/,  is  only  first-order  accurate. 

Assume  that  at  any  time  level  n  the  divergence-free  body  velocity  vb  is  known.  For  body  cells 


with  (Jv  —  7  ,  the  no-slip  condition  requires  that  v  =  vb  .  For  partially  filled  interface  cells  with 


0  <§b  <1 ,  special  care  must  be  taken  to  account  for  the  influences  of  the  body  surface  as  accurately 


as  possible.  In  Xiao  [10],  the  cell-center  velocity  is  defined  as  the  volume  average  of  the  body 
velocity  and  the  fluid  velocity  computed  for  the  particular  interface  cell.  In  Fadlun  et  al.  [2]  and  Kim 
et  al.  [3],  the  velocity  of  the  interface  cells  is  not  computed  by  the  flow  equations  but  interpolated 
using  the  surrounding  fluid  and  body  cells.  In  our  experience,  both  approaches  will  work  when  the 
immersed  body  is  stationary.  However,  when  the  body  is  accelerating,  an  impulsive  change  in 
pressure  will  be  created  corresponding  to  the  change  in  body  velocity.  Which  approach  can  model  this 
pressure  change  correctly  requires  further  investigation  efforts.  Here,  we  show  some  results  using  the 
volume-averaging  approach. 

The  direct  momentum  forcing  in  Eq.  (2)  is  written  as 


A  -  -§t 


f  —  —  n 

vb-v 


V 


At 


AV+9(RS  +  J^R'"AS)  +  (l-9)(Rn  +^P"AX) 


\ 


(15) 


cs  cs 

Thus,  for  fluid  cells  with  <§>h  =  0 ,  no  forcing  is  applied.  For  body  cells  with  (\>b  =  7 ,  the  converged 


_ k  i):  _ k 

solution  of  Eq.  (2)  is  v  =  vh  . 


For  interface  cells  with  0  <(. \>h  <  1 ,  the  solution  is  a  volume  average 


of  the  body  velocity  vb  and  the  fluid  velocity  v  ,Iuid  .  However,  the  fluid  velocity  v  fhiid  is  not 

explicitly  computed.  Instead,  the  volume-averaged  v  computed  by  Eq.  (2)  is  used  as  the  fluid 
velocity. 

As  for  the  corrections  for  velocity  and  pressure,  we  assume  that  Eqs.  (5),  (6),  (8)  and  (9)  are  valid 
for  interface  cells.  For  body  cells,  there  is  no  need  for  velocity  correction.  Furthermore,  because  the 
body  velocity  is  assumed  divergence  free,  it  is  reasonable  to  assume  that  inside  the  body  cells  the 
pressure  correction  satisfies  Eq.  (6)  with  a  zero  source  term,  or  equivalently  the  Laplace  equation.  In 
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essence  this  is  equivalent  to  treating  the  body  as  an  incompressible  fluid  with  prescribed  velocity 
distribution,  and  hence  it  satisfies  the  same  pressure  correction  equation  as  the  fluid  outside  the  body. 

9.  Driven  Cavity  With  Centered  Cylinders 

A  stationary  circular  cylinder  of  radius  0.21  is  placed  at  the  center  of  the  driven  cavity  in  the  last 
example.  The  cylinder  radius  is  arbitrarily  chosen.  The  grid  has  one  level  of  refinement  around  the 
cylinder,  as  shown  in  Fig.  6.  The  Reynolds  number  based  on  the  wall  speed  and  the  cavity  width  is 
1000.  The  computed  particle  paths  are  plotted  in  Fig.  7.  The  contour  of  VOB  function  (f)&  = 0.5  is 
also  plotted  as  the  cylinder  wall.  The  smoothness  of  the  streamlines  over  the  cylinder  wall  indicates 
that  the  treatment  of  the  interface  cells  is  effective  in  representing  a  curved  surface.  For  comparison’ s 
sake,  a  multi-zone  body-fitted  flow  solver  [6]  is  used  to  solve  the  same  problem  on  a  5-zone  grid.  The 
velocity  component  u(y)  along  the  vertical  centerline  and  v(x)  along  horizontal  centerline  are 
compared  in  Fig.  8.  The  lines  are  computed  by  the  body-fitted  method  and  the  symbols  are  computed 
by  the  immersed  boundary  method.  Because  the  Although  some  difference  can  be  seen  in  the  figure, 
but  generally  the  shapes  of  the  velocity  profiles  compare  very  well. 

The  pressure  and  viscous  forced  experienced  by  the  cylinder  are  also  compared.  For  the  body- 
fitted  method,  the  pressure  and  viscous  stresses  are  integrated  along  the  cylinder  surface  to  obtain  the 
forces.  For  the  current  method,  the  forces  experienced  by  the  body  are  obtained  by 

A*  =  E4>„(-VP  +  2-V2v)AV  (16) 

all  & e 

The  body-fitted  method  obtained  a  pressure  force  of  (-0.4619E-2,  -0.4418E-4)  and  a  viscous  force  of 
(-0.1734E-2,  0.6231E-3).  The  current  method  obtained  a  pressure  force  of  (-0.4619E-2,  -0.1995E-3) 
and  a  viscous  force  of  (-0.1305E-2,  0.4773E-3).  The  pressure  forces  are  in  better  agreement  than  the 
viscous  forces  are.  This  is  not  surprising,  since  the  immersed  boundary  method  mainly  influences  the 
velocity  profde  near  body  surfaces  and  the  pressure  force  is  mainly  an  inviscid  effect. 

10.  Conclusions 

An  implicit  pressure  correction  method  for  the  incompressible  Navier-Stokes  equations  on 
unstructured  Cartesian  grid  is  developed  and  tested.  The  technique  of  sub -iteration  is  used  to  derive 
the  implicit  time  integration  equation.  The  pressure  Poisson  equation  is  solved  by  an  implicit  multi¬ 
grid  method.  A  forth-order  dissipation  is  added  to  the  pressure  field  to  control  the  even-odd 
decoupling.  Tests  show  that  with  an  appropriate  amount  of  artificial  dissipation,  the  present  method  is 
second-order  accurate  both  in  time  and  space.  An  immersed  boundary  method  is  also  developed  and 
tested.  The  results  indicate  that  the  current  method  is  simple  and  effective  in  modeling  the  stationary 
immersed  bodies.  The  extended  use  of  the  present  approach  to  moving  bodies  is  currently  undergoing. 
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dx  dt 


Fig.  1  Spatial  accuracy  test,  decaying  vortices,  Fig.  3  Time  accuracy  test,  decaying  vortices, 
cdisp=0.03  0  =0.5 


Fig.  2  Spatial  accuracy  test,  decaying  vortices,  Fig.  4  Driven  Cavity,  u(y)  along  x=0.5  and  v(x) 
cdisp=0.  along  y=0.5,  Re=3200,  cdisp=0.01 
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X 

Fig.  8  Velocity  u(y)  along  x=0.5  and  v(x)  along 
Fig.  5  Particle  paths,  Re=3200  v=0.5,  driven  cavity  with  centered  cylinder 


0.25  0.5  0.75 

X 

Fig.  6  Unstructured  Cartesian  grid  for  cavity  with 
embedded  cylinder,  one  level  refinement 


Fig.  7  Particle  Paths  for  driven  cavity  with 
centered  cylinder,  Re=1000 
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